Analisar e comparar o desempenho de diferentes técnicas de modelagem de séries temporais na previsão da Série de Empregos Formais em Minas Gerais, visando atender os 4 pontos abaixo.
Em Anexo, estará todo o codigo utilizado
A Regressão em Séries Temporais é uma técnica estatística que estabelece uma relação linear entre uma variável dependente (a série temporal a ser prevista) e uma ou mais variáveis independentes (fatores explicativos). Através de um modelo de regressão linear, a técnica busca identificar como as variáveis independentes influenciam o comportamento da série temporal ao longo do tempo. A Regressão em Séries Temporais é uma técnica para a análise e previsão de séries temporais, mas deve ser utilizada com cautela, considerando suas suposições e limitações.
Após 5 tentativas anteriores, o modelo escolhido foi a sexta versao, visto que apresentou um menor AIC, o melhor r quadrado ajustado, e menor termo de erro. Este modelo é composto pelo tempo, pelo quadrado do tempo (para explicar a tendencia quadrática observada na serie), pelo fator sazonal e por um termo autoregressivo de ordem 1. Com isto, tem-se os pontos abaixo:
Modelo 6 - Termo autoregressivo na serie transformada
resid_1=rep(0,nEmp)
for(i in 2:nEmp)
resid_1[i]=modelo5$residuals[i-1]
modelo6= lm(lnEmp ~ t + t2 + factor(sazon) + resid_1)
summary(modelo6)
##
## Call:
## lm(formula = lnEmp ~ t + t2 + factor(sazon) + resid_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.005222 -0.001443 0.000167 0.001492 0.004962
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.504e+00 1.537e-03 2930.595 < 2e-16 ***
## t -2.929e-04 7.098e-05 -4.127 0.000139 ***
## t2 3.774e-05 1.046e-06 36.071 < 2e-16 ***
## factor(sazon)2 2.089e-03 1.631e-03 1.281 0.206195
## factor(sazon)3 2.176e-03 1.631e-03 1.334 0.188363
## factor(sazon)4 8.466e-03 1.569e-03 5.397 1.86e-06 ***
## factor(sazon)5 2.190e-02 1.568e-03 13.967 < 2e-16 ***
## factor(sazon)6 3.690e-02 1.568e-03 23.540 < 2e-16 ***
## factor(sazon)7 3.945e-02 1.568e-03 25.164 < 2e-16 ***
## factor(sazon)8 3.596e-02 1.568e-03 22.933 < 2e-16 ***
## factor(sazon)9 3.002e-02 1.645e-03 18.252 < 2e-16 ***
## factor(sazon)10 2.046e-02 1.632e-03 12.536 < 2e-16 ***
## factor(sazon)11 1.624e-02 1.631e-03 9.956 1.87e-13 ***
## factor(sazon)12 1.444e-03 1.631e-03 0.885 0.380261
## resid_1 9.718e-01 7.136e-02 13.619 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002579 on 50 degrees of freedom
## Multiple R-squared: 0.9975, Adjusted R-squared: 0.9968
## F-statistic: 1434 on 14 and 50 DF, p-value: < 2.2e-16
AIC(modelo6)
## [1] -575.4465
Teste de normalidade dos residuos
plot(modelo6)
shapiro.test(modelo6$residuals)
##
## Shapiro-Wilk normality test
##
## data: modelo6$residuals
## W = 0.98736, p-value = 0.7498
Não ha evidencias para rejeitar H0 (H0: os residuos sao normalmente distribuidos)
durbinWatsonTest(modelo6)
## lag Autocorrelation D-W Statistic p-value
## 1 0.3358639 1.279857 0.004
## Alternative hypothesis: rho != 0
Há evidencias para rejeitar-se H0 - H0: nao ha auto-correlacao de ordem 1
acf(modelo6$residuals)
Ratifica-se o valor observado no teste de Durbin-Watson - ve-se
autocorrelações significante
Foi removido os ultimos 9 registros da série para fazer a previsão e comparação entre previsto e realizado.
real=lnEmp[57:65]
n6=length(lnEmp)-9 # Coloca em n o tamanho da serie menos 12 observacoes
serie=rep(0,n6) # Cria a variavel serie
for(i in 1:n6)
serie[i]=lnEmp[i] # Coloca na variavel serie a serie transformada por LN
serie=ts(serie,start=c(1999,04),frequency=12)
t=c(1:n6)
t2=t^2
sazon1=rep.int(4:12,1)
sazon2=rep.int(1:12,3)
sazon3=rep.int(1:11,1)
sazon6 = append(sazon1, sazon2)
sazon6 = append(sazon6, sazon3)
M5=lm(serie ~ t + t2+ factor(sazon6))
resid_m5=rep(0,n6)
for(i in 2:n6)
resid_m5[i]=M5$res[i-1]
M6= lm(serie ~ t + t2 + factor(sazon6) + resid_m5)
summary(M6)
##
## Call:
## lm(formula = serie ~ t + t2 + factor(sazon6) + resid_m5)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0048750 -0.0012220 -0.0001147 0.0011958 0.0052264
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.502e+00 1.629e-03 2763.130 < 2e-16 ***
## t 3.143e-05 8.442e-05 0.372 0.712
## t2 3.148e-05 1.437e-06 21.909 < 2e-16 ***
## factor(sazon6)2 2.185e-03 1.734e-03 1.260 0.215
## factor(sazon6)3 2.153e-03 1.734e-03 1.242 0.221
## factor(sazon6)4 8.771e-03 1.655e-03 5.300 4.25e-06 ***
## factor(sazon6)5 2.151e-02 1.653e-03 13.013 3.80e-16 ***
## factor(sazon6)6 3.607e-02 1.652e-03 21.829 < 2e-16 ***
## factor(sazon6)7 3.763e-02 1.652e-03 22.779 < 2e-16 ***
## factor(sazon6)8 3.310e-02 1.652e-03 20.038 < 2e-16 ***
## factor(sazon6)9 2.695e-02 1.653e-03 16.310 < 2e-16 ***
## factor(sazon6)10 2.013e-02 1.654e-03 12.173 3.38e-15 ***
## factor(sazon6)11 1.599e-02 1.656e-03 9.656 4.06e-12 ***
## factor(sazon6)12 1.968e-03 1.734e-03 1.135 0.263
## resid_m5 8.793e-01 7.456e-02 11.794 9.34e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002451 on 41 degrees of freedom
## Multiple R-squared: 0.9962, Adjusted R-squared: 0.9949
## F-statistic: 765.8 on 14 and 41 DF, p-value: < 2.2e-16
previsao = rep(1,9)
previsao[1] = M6$coef[1] + M6$coef[2]*(n6+1) + M6$coef[3]*(n6+1)^2 + M6$coef[14] + M6$coef[15]*M5$res[n6]
previsao[2] = M6$coef[1] + M6$coef[2]*(n6+2) + M6$coef[3]*(n6+2)^2 + M6$coef[15]*M5$res[n6]
for(i in 3:9)
previsao[i] = M6$coef[1] + M6$coef[2]*(n6+i) + M6$coef[3]*(n6+i)^2 + M6$coef[i+1] + (M6$coef[15])^i*M6$res[n6]
previsao
## [1] 4.608653 4.610338 4.615821 4.619523 4.629942 4.646553 4.665044 4.670609
## [9] 4.670153
exp(previsao)
## [1] 100.3489 100.5181 101.0708 101.4456 102.5082 104.2252 106.1703 106.7627
## [9] 106.7141
previsao=ts(previsao,start = c(2003,12), frequency = 12)
plot(exp(previsao))
real=ts(real,start = c(2003,12), frequency = 12)
plot(exp(real))
previsao
## Jan Feb Mar Apr May Jun Jul Aug
## 2003
## 2004 4.610338 4.615821 4.619523 4.629942 4.646553 4.665044 4.670609 4.670153
## Sep Oct Nov Dec
## 2003 4.608653
## 2004
real
## Jan Feb Mar Apr May Jun Jul Aug
## 2003
## 2004 4.609162 4.615121 4.620059 4.630838 4.652054 4.673763 4.685828 4.692265
## Sep Oct Nov Dec
## 2003 4.605170
## 2004
Intervalos de previsao
li=rep(1,9)
ls=rep(1,9)
previsaofinal=rep(1,9)
p=length(M6$coef)
PHI=1
for(i in 1:9)
{
li[i]=previsao[i] - qt(0.975,n6-p)*sd(M6$res)*sqrt(PHI)
ls[i]=previsao[i] + qt(0.975,n6-p)*sd(M6$res)*sqrt(PHI)
PHI=PHI+(M6$coef[15])^2*i
}
li
## [1] 4.604379 4.604646 4.608033 4.609372 4.617311 4.631382 4.647297 4.650265
## [9] 4.647200
for(i in 1:9)
{
previsaofinal[i] = exp(previsao[i])
li[i]=exp(li[i])
ls[i]=exp(ls[i])
}
length(previsaofinal)
## [1] 9
length(real)
## [1] 9
vetor_prev = cbind(exp(real),previsaofinal, li, ls)
vetor_prev
## exp(real) previsaofinal li ls
## Dec 2003 100.0 100.3489 99.92087 100.7788
## Jan 2004 100.4 100.5181 99.94755 101.0919
## Feb 2004 101.0 101.0708 100.28673 101.8610
## Mar 2004 101.5 101.4456 100.42111 102.4806
## Apr 2004 102.6 102.5082 101.22151 103.8112
## May 2004 104.8 104.2252 102.65580 105.8185
## Jun 2004 107.1 106.1703 104.30270 108.0713
## Jul 2004 108.4 106.7627 104.61275 108.9568
## Aug 2004 109.1 106.7141 104.29251 109.1918
Calculo do Erro Quadratico Medio de Previsao:
for(i in 1:9)
SQP=(exp(real[i]) - previsaofinal[i])^2/9
SQP
## [1] 0.6325209
O valor observado acima é pequeno, o que mostra uma aderencia boa para o modelo de regressão em séries temporais.
O modelo SARIMA (Seasonal Autoregressive Integrated Moving Average) combina a robustez do modelo ARIMA com a flexibilidade da Suavização Exponencial para prever séries temporais com sazonalidade. Essa técnica poderosa se destaca pela capacidade de capturar padrões sazonais complexos e adaptar-se a mudanças no comportamento da série ao longo do tempo. O modelo SARIMA é utilizado para previsão de séries temporais com sazonalidade. Sua capacidade de capturar padrões sazonais complexos e se adaptar a mudanças no comportamento da série o torna uma escolha popular para diversas aplicações. No entanto, é importante considerar as suposições do modelo e a necessidade de uma seleção cuidadosa dos parâmetros para garantir a precisão das previsõ
Foram feitos alguns testes de melhor ajuste dos parametros Auto Regressivo, Diferencial, Média Movel e Sazonal. Após as analises das saidas do modelo, observa-se que o modelo de melhor ajuste para a serie do Emprego Formal de Minas Gerais, foi o modelo abaixo, que apresentou o menor AIC.
fit2 <- sarima(lnEmp,1,1,0,1,1,1,6)
## initial value -4.129883
## iter 2 value -5.469381
## iter 3 value -5.566620
## iter 4 value -5.571658
## iter 5 value -5.576403
## iter 6 value -5.578365
## iter 7 value -5.578528
## iter 8 value -5.578827
## iter 9 value -5.578846
## iter 10 value -5.578851
## iter 11 value -5.578854
## iter 12 value -5.578859
## iter 12 value -5.578859
## iter 12 value -5.578859
## final value -5.578859
## converged
## initial value -5.475502
## iter 2 value -5.478572
## iter 3 value -5.482324
## iter 4 value -5.488341
## iter 5 value -5.489446
## iter 6 value -5.489628
## iter 7 value -5.489638
## iter 8 value -5.489640
## iter 8 value -5.489640
## final value -5.489640
## converged
## <><><><><><><><><><><><><><>
##
## Coefficients:
## Estimate SE t.value p.value
## ar1 0.5956 0.1249 4.7683 0.0000
## sar1 -0.8863 0.0539 -16.4451 0.0000
## sma1 -0.5848 0.2012 -2.9066 0.0053
##
## sigma^2 estimated as 1.27053e-05 on 55 degrees of freedom
##
## AIC = -8.003472 AICc = -7.995809 BIC = -7.861373
##
fit2
## $fit
##
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S),
## include.mean = !no.constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 sar1 sma1
## 0.5956 -0.8863 -0.5848
## s.e. 0.1249 0.0539 0.2012
##
## sigma^2 estimated as 1.271e-05: log likelihood = 236.1, aic = -464.2
##
## $degrees_of_freedom
## [1] 55
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.5956 0.1249 4.7683 0.0000
## sar1 -0.8863 0.0539 -16.4451 0.0000
## sma1 -0.5848 0.2012 -2.9066 0.0053
##
## $ICs
## AIC AICc BIC
## -8.003472 -7.995809 -7.861373
shapiro.test(fit2$fit$residuals)
##
## Shapiro-Wilk normality test
##
## data: fit2$fit$residuals
## W = 0.96492, p-value = 0.06214
Nota-se que o modelo não possui media movel para componente nao sazonal, e isto devido ao fato desta componente não apresentar significancia na analise da versão anterior do modelo. Através da análise da saída do modelo, é possivel comprovar a condição de inversibilidade, para o componente MA sazonal, e também é possivel comprovar a estacionariedade dos componentes autoregressivos, seja sazonal ou não sazonal.
Da mesma maneira que a técnica anterior, exclui-se as 9 ultimas observações e faz-se a comparação entre o real e o previsto, para obtenção do indicador de qualidade de ajuste (MSE). Importante ressaltar que a série utilizada para o ajuste do SARIMA também está transformada por logaritmo.
serie
## Jan Feb Mar Apr May Jun Jul Aug
## 1999 4.511958 4.526127 4.539030 4.541165 4.536891
## 2000 4.503137 4.506454 4.506454 4.514151 4.528289 4.548600 4.555980 4.554929
## 2001 4.518522 4.521789 4.522875 4.530447 4.541165 4.555980 4.558079 4.553877
## 2002 4.535820 4.541165 4.546481 4.557030 4.575741 4.591071 4.594109 4.593098
## 2003 4.574711 4.578826 4.579852 4.587006 4.601162 4.619073 4.621044 4.618086
## Sep Oct Nov Dec
## 1999 4.526127 4.519612 4.517431 4.503137
## 2000 4.549657 4.537961 4.534748 4.520701
## 2001 4.549657 4.548600 4.544358 4.533674
## 2002 4.591071 4.585987 4.584967 4.574711
## 2003 4.619073 4.619073 4.619073
serieSarima <- sarima.for(serie, 9, 1, 1, 0, 1, 1, 1, 6)
serieSarima$pred
## Jan Feb Mar Apr May Jun Jul Aug
## 2003
## 2004 4.609670 4.614122 4.613904 4.618848 4.631834 4.647129 4.649037 4.646800
## Sep Oct Nov Dec
## 2003 4.608695
## 2004
exp(serieSarima$pred)
## Jan Feb Mar Apr May Jun Jul Aug
## 2003
## 2004 100.4510 100.8992 100.8772 101.3772 102.7023 104.2851 104.4843 104.2508
## Sep Oct Nov Dec
## 2003 100.3531
## 2004
Novamente calcula-se o erro quadrático médio
for(i in 1:9)
SQP=(exp(real[i]) - exp(serieSarima$pred[i]))^2/9
SQP
## [1] 2.61274
E já é possivel notar que embora seja usada uma técnica mais robusta, o ajuste deste modelo SARIMA ficou pior que o ajuste do modelo de regressão em séries temporais, observando o MSE de ambos.
A Suavização Exponencial é uma técnica de previsão robusta e adaptável que se destaca pela simplicidade e eficiência. Através da atribuição de pesos exponenciais decrescentes às observações passadas, a técnica captura as tendências recentes da série temporal e gera previsões precisas, especialmente para séries com ruído e mudanças bruscas
De fato a implementação computacional em código R para esta técnica é muito simples, bastando apenas o uso de uma função.
fitHT <- HoltWinters(lnEmp)
plot(fitted(fitHT))
plot(fitHT)
Observando o gráfico gerado pela saída do modelo, é possivel observar os 3 componentes, um nivel, uma tendencia e uma sazonalidade, além de confirmar o bom ajuste, linha vermelha, seguindo a linha preta.
Seguindo o padrão ja realizado anteriormente neste relatório, utiliza-se o vetor com as ultimas 9 observações removidas, ajusta-se um novo modelo para fazer a comparação com o realizado.
fitHT2 <- HoltWinters(serie)
plot(fitted(fitHT2))
plot(fitHT2)
predHT2 = predict(fitHT2,9, prediction.interval = TRUE)
exp(predHT2)
## fit upr lwr
## Dec 2003 99.83901 100.6652 99.01956
## Jan 2004 99.77359 100.7571 98.79969
## Feb 2004 100.18276 101.3186 99.05965
## Mar 2004 100.38318 101.6638 99.11871
## Apr 2004 101.26748 102.6992 99.85568
## May 2004 102.81448 104.4077 101.24562
## Jun 2004 104.60144 106.3628 102.86927
## Jul 2004 105.04856 106.9577 103.17351
## Aug 2004 105.10032 107.1504 103.08943
exp(real)
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2003 100.0
## 2004 100.4 101.0 101.5 102.6 104.8 107.1 108.4 109.1
plot(exp(predHT2))
plot(exp(real))
Com o modelo ajustado, é possivel avaliar o erro médio quadratico.
for(i in 1:9)
SQP=(exp(real[i]) - exp(predHT2[i,1]))^2/9
SQP
## fit
## 1.77749
E o valor observado foi menor que o da técnica SARIMA, porém ainda é maior que o erro médio quadratico encontrado usando a técnica de regressão em séries temporais.
O modelo de Espaço-Estado se destaca por sua flexibilidade e capacidade de representar sistemas dinâmicos complexos com alta precisão. Através de um sistema de equações diferenciais e matrizes, a técnica modela as relações entre variáveis observáveis e não observáveis do sistema, capturando a dinâmica interna e as interações entre os componentes. Essa característica o torna ideal para prever o comportamento futuro de sistemas que apresentam mudanças não lineares, sazonalidade e ruído. é importante considerar a complexidade da implementação, a necessidade de expertise e a disponibilidade de dados de alta qualidade para garantir a efetividade da técnica. É importante considerar a complexidade da implementação, a necessidade de expertise e a disponibilidade de dados de alta qualidade para garantir a efetividade da técnica.
Como colocado na sessão anterior, é complexo modelar usando sistemas de espaço-estado, e o código em R também é mais complexo. É necessário ter maior dominio e maior desenvolvimento de código. E para esta técnica, a série transformada gera erro na saida do modelo, ou seja, foi necessário utilizar a série original.
model<-function(u){
mod<-dlmModSeas(frequency=12,dV=0,dW=c(exp(u[4]),rep(0,10)))+
dlmModPoly(2,dV=exp(u[1]),dW=(exp(u[2:3])))
}
outmle=dlmMLE(tsEmp,parm=rep(0,4),model)
exp(outmle$par)
## [1] 1.667719e-06 3.801566e-02 1.067618e-02 3.517537e-08
mod=model(outmle$par)
outmodFil=dlmFilter(tsEmp,mod)
outF<-dlmFilter(tsEmp,mod) #Filtering
outS<-dlmSmooth(tsEmp,mod) #smoothing
par(mfrow=c(4,1))
ts.plot((outF$m[10:nEmp,1]))
title("Sazonality: Filtering")
ts.plot((outF$m[10:nEmp,12]))
title("Slope: Filtering")
ts.plot((outF$m[10:nEmp,13]))
title("Level: Filtering")
ts.plot((outF$f[10:nEmp]),ylab="Yt^",xlab="t")
title("Preditos")
par(mfrow=c(1,1))
qqnorm(residuals(outF,sd=FALSE))
qqline(residuals(outF,sd=FALSE))
tsdiag(outF)
shapiro.test(residuals(outF,sd=FALSE))
##
## Shapiro-Wilk normality test
##
## data: residuals(outF, sd = FALSE)
## W = 0.94949, p-value = 0.009997
Tendo em vista os 4 graficos gerados pelo codigo acima, nota-se que os primeiros pontos para todos eles, estão, de certa forma, um pouco fora da escala, enquanto que os demais, vão, de certa forma mais alinhados. Uma caracteristica de adequação de sistemas dinâmicos na engenharia. Importante registar que o teste de Sharipo-wilk feito em cima dos residuos do modelo de Espaço-Estado, tem um p-valor que traz evidências para rejeitar H0, ou seja, os residuos não são normalmente distribuidos.
Mesmo com a ressalva acima, dos residuos do modelo não ser normalmente distribuido, foi definido fazer a previsão usando o modelo ajustado para a série sem as ultimas 9 observações.
seriePura = exp(serie)
seriePura
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1999 91.1 92.4 93.6 93.8 93.4 92.4 91.8 91.6 90.3
## 2000 90.3 90.6 90.6 91.3 92.6 94.5 95.2 95.1 94.6 93.5 93.2 91.9
## 2001 91.7 92.0 92.1 92.8 93.8 95.2 95.4 95.0 94.6 94.5 94.1 93.1
## 2002 93.3 93.8 94.3 95.3 97.1 98.6 98.9 98.8 98.6 98.1 98.0 97.0
## 2003 97.0 97.4 97.5 98.2 99.6 101.4 101.6 101.3 101.4 101.4 101.4
outmle2=dlmMLE(seriePura,parm=rep(0,4),model)
exp(outmle2$par)
## [1] 1.948774e-06 4.435990e-02 3.268541e-03 1.822108e-07
mod2=model(outmle2$par)
outmodFil2=dlmFilter(seriePura,mod2)
outF2<-dlmFilter(seriePura,mod2) #Filtering
outS2<-dlmSmooth(seriePura,mod2) #smoothing
prev2=dlmForecast(outmodFil2,n=9)
prev2$f
## Jan Feb Mar Apr May Jun Jul Aug
## 2003
## 2004 100.7009 101.2891 101.6691 102.6409 104.2304 106.0133 106.5496 106.4995
## Sep Oct Nov Dec
## 2003 100.4796
## 2004
previsaoMEB=as.numeric(prev2$f)
Execução do MSE, para comparação com os demais modelos gerados anteriormente.
for(i in 1:9)
SQP=(exp(real[i]) - previsaoMEB[i])^2/9
SQP
## [1] 0.7514231
O valor do MSE observado pelo modelo ajustado usando a técnica de Espaço-Estado foi de 0,75, o segundo melhor gerado neste estudo.
Para a série temporal de empregos formais no estado de Minas Gerais entre abril de 1999 e agosto de 2004, conclui-se que o modelo usando a técnica de regressão em séries temporais, se mostrou o mais precisa. Embora o esforço para se chegar ao modelo adequado, através desta técnica, tenha sido o maior, como pode ser visto no Anexo. Em seguida, um modelo de complexa implementação de Espaço-Estado, teve o segundo melhor desepenho em precisão, seguido depois pelo modelo de suavização exponencial, o mais simples de implementar técnicamente, e por fim, o modelo SARIMA, teve o pior desempenho para prever as ultimas 9 observações com precisão. Porém, em todos os casos os desvios foram pequenos, desta forma, para está série, pensando em esforço e retorno, a técnica de suavização exponencial seria a mais adequada a ser implementada, dada a simplicidade, sem perder precisão de forma significativa.
library(FNN)
library(corrplot)
library(ggplot2)
require(car)
library(tidyverse)
library(ks)
library(tseries)
library(moments)
library(stats)
library(forecast)
library(astsa)
library(dlm)
setwd("/Users/victortelles/Documents/Coursera/Especializacao - UFMG/06 - Analise de Series Temporais/Trabalho Pratico Final/serie2")
stEmp <- read.csv("SerieIEF.csv", sep=";", header = T)
nEmp = length(stEmp$yt)
ts.plot(stEmp$yt, xlab = "tempo", ylab = "Indice de emprego formal em MG")
tsEmp <- ts(stEmp$yt, start = c(1999,4), frequency = 12)
tsEmp
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1999 91.1 92.4 93.6 93.8 93.4 92.4 91.8 91.6 90.3
## 2000 90.3 90.6 90.6 91.3 92.6 94.5 95.2 95.1 94.6 93.5 93.2 91.9
## 2001 91.7 92.0 92.1 92.8 93.8 95.2 95.4 95.0 94.6 94.5 94.1 93.1
## 2002 93.3 93.8 94.3 95.3 97.1 98.6 98.9 98.8 98.6 98.1 98.0 97.0
## 2003 97.0 97.4 97.5 98.2 99.6 101.4 101.6 101.3 101.4 101.4 101.4 100.0
## 2004 100.4 101.0 101.5 102.6 104.8 107.1 108.4 109.1
plot(tsEmp)
## Analise Exploratoria da Série Temporal
cv = sd(tsEmp)/mean(tsEmp)
cv
## [1] 0.04655599
Em media a série varia em 4,7%
skewness(tsEmp)
## [1] 0.8733008
kurtosis(tsEmp)
## [1] 3.275672
Observa-se simetria em torno da média. Pela curtose, temos possibilidade de gerar bastante valores nas caldas para valores de Curtoses alto - valor de referencia 3
summary(tsEmp)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 90.30 93.10 95.10 96.35 98.90 109.10
boxplot(tsEmp)
fatorSazonal=rep(4:12,1)
fatorSazonal2=rep(1:12,4)
fatorSazonal3=rep(1:8,1)
fatorSazo = append(fatorSazonal, fatorSazonal2)
fatorSazo = append(fatorSazo, fatorSazonal3)
boxplot(tsEmp~fatorSazo,xlab="Meses",ylab="Ãndice de emprego formal em MG")
Observando o boxplot mensal, é possivel observar as diferenças
pluviometricas ao longo do ano
hist(tsEmp)
Ratifica-se o valor de curtose, observando caldas longas
qqnorm(tsEmp); qqline(tsEmp)
Observando o qqplot, pode-se dizer que os dados não respeitam a
distribuição normal
shapiro.test(tsEmp)
##
## Shapiro-Wilk normality test
##
## data: tsEmp
## W = 0.92401, p-value = 0.0006615
Há evidencias para rejeitar H0 (H0: de que os dados sao normalmente distribuidos)
acf(tsEmp)
Observa-se que a serie tem memoria longa - alem de vermos as ondas
tradicionais de sazonalidade
pacf(tsEmp)
Observa-se que a ordem 1 ultrapassa as bandas limites
Box.test(tsEmp, lag = 1, type="Ljung-Box")
##
## Box-Ljung test
##
## data: tsEmp
## X-squared = 55.865, df = 1, p-value = 7.76e-14
Há evidencias para rejeitar H0 (H0: o conjunto de correlação de ordem 1 é igual zero)
Box.test(tsEmp, lag = 6, type="Ljung-Box")
##
## Box-Ljung test
##
## data: tsEmp
## X-squared = 201.28, df = 6, p-value < 2.2e-16
Há evidencias para rejeitar H0 (H0: o conjunto de correlação de ordem 6 é igual zero) ou seja possui memória de mais de 6 meses
Box.test(tsEmp, lag = 12, type="Ljung-Box")
##
## Box-Ljung test
##
## data: tsEmp
## X-squared = 309.24, df = 12, p-value < 2.2e-16
Há evidencias para rejeitar H0 (H0: o conjunto de correlação de ordem 12 é igual zero) ou seja possui memória de mais de 12 meses
Box.test(tsEmp, lag = 24, type="Ljung-Box")
##
## Box-Ljung test
##
## data: tsEmp
## X-squared = 344.87, df = 24, p-value < 2.2e-16
Há evidencias para rejeitar H0 (H0: o conjunto de correlação de ordem 24 é igual zero) ou seja possui memória de mais de 24 meses
Box.test(tsEmp, lag = 36, type="Ljung-Box")
##
## Box-Ljung test
##
## data: tsEmp
## X-squared = 411.54, df = 36, p-value < 2.2e-16
Há evidencias para rejeitar H0 (H0: o conjunto de correlação de ordem 36 é igual zero) ou seja possui memória de mais de 36 meses
t = c(1:nEmp)
modelo1 <- lm(tsEmp ~ t)
summary(modelo1)
##
## Call:
## lm(formula = tsEmp ~ t)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2653 -1.4309 -0.3491 1.0658 5.9786
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 89.37481 0.51424 173.80 <2e-16 ***
## t 0.21149 0.01355 15.61 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.049 on 63 degrees of freedom
## Multiple R-squared: 0.7946, Adjusted R-squared: 0.7913
## F-statistic: 243.7 on 1 and 63 DF, p-value: < 2.2e-16
AIC(modelo1)
## [1] 281.6922
plot(modelo1)
shapiro.test(modelo1$residuals)
##
## Shapiro-Wilk normality test
##
## data: modelo1$residuals
## W = 0.95311, p-value = 0.01517
Há evidencias para rejeitar H0 (H0: os residuos sao normalmente distribuidos)
durbinWatsonTest(modelo1)
## lag Autocorrelation D-W Statistic p-value
## 1 0.8390991 0.1780147 0
## Alternative hypothesis: rho != 0
Há evidencias para rejeitar-se H0 - H0: nao ha auto-correlacao de ordem 1
acf(modelo1$residuals)
Ratifica-se o valor observado no teste de Durbin-Watson - vê-se
autocorrelações significantes
t2=t^2
modelo2 <- lm(tsEmp ~ t + t2)
summary(modelo2)
##
## Call:
## lm(formula = tsEmp ~ t + t2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.42808 -1.23632 -0.04589 1.18878 3.05655
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 92.5795009 0.5825264 158.928 < 2e-16 ***
## t -0.0755014 0.0407283 -1.854 0.0685 .
## t2 0.0043483 0.0005981 7.270 7.31e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.518 on 62 degrees of freedom
## Multiple R-squared: 0.8891, Adjusted R-squared: 0.8856
## F-statistic: 248.6 on 2 and 62 DF, p-value: < 2.2e-16
AIC(modelo2)
## [1] 243.615
plot(modelo2)
shapiro.test(modelo2$residuals)
##
## Shapiro-Wilk normality test
##
## data: modelo2$residuals
## W = 0.96047, p-value = 0.03629
Há evidencias para rejeitar-se H0 - H0: nao ha auto-correlacao de ordem 1
acf(modelo2$residuals)
Ratifica-se o valor observado no teste de Durbin-Watson - vê-se autocorrelações significantes
sazon1=rep.int(4:12,1)
sazon2=rep.int(1:12,4)
sazon3=rep.int(1:8,1)
sazon = append(sazon1, sazon2)
sazon = append(sazon, sazon3)
modelo3 <- lm(tsEmp ~ t + t2 + factor(sazon))
summary(modelo3)
##
## Call:
## lm(formula = tsEmp ~ t + t2 + factor(sazon))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.21429 -0.36498 -0.04533 0.33957 1.69858
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 90.3224275 0.3473353 260.044 < 2e-16 ***
## t -0.0352069 0.0159769 -2.204 0.032093 *
## t2 0.0037497 0.0002349 15.966 < 2e-16 ***
## factor(sazon)2 0.1964757 0.3688266 0.533 0.596552
## factor(sazon)3 0.2054519 0.3688916 0.557 0.580002
## factor(sazon)4 0.8072789 0.3547153 2.276 0.027086 *
## factor(sazon)5 2.1062529 0.3545473 5.941 2.54e-07 ***
## factor(sazon)6 3.5810608 0.3544634 10.103 9.09e-14 ***
## factor(sazon)7 3.8483693 0.3544627 10.857 7.26e-15 ***
## factor(sazon)8 3.5248450 0.3545462 9.942 1.57e-13 ***
## factor(sazon)9 2.5991027 0.3691354 7.041 4.69e-09 ***
## factor(sazon)10 1.9455762 0.3689891 5.273 2.75e-06 ***
## factor(sazon)11 1.5445503 0.3688869 4.187 0.000112 ***
## factor(sazon)12 0.1360249 0.3688260 0.369 0.713800
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5831 on 51 degrees of freedom
## Multiple R-squared: 0.9865, Adjusted R-squared: 0.9831
## F-statistic: 287.4 on 13 and 51 DF, p-value: < 2.2e-16
AIC(modelo3)
## [1] 128.5813
plot(modelo3)
shapiro.test(modelo3$residuals)
##
## Shapiro-Wilk normality test
##
## data: modelo3$residuals
## W = 0.97801, p-value = 0.2995
Não ha evidencias para rejeitar H0 (H0: os residuos sao normalmente distribuidos)
durbinWatsonTest(modelo3)
## lag Autocorrelation D-W Statistic p-value
## 1 0.7918974 0.249839 0
## Alternative hypothesis: rho != 0
há evidencias para rejeitar-se H0 - H0: nao ha auto-correlacao de ordem 1
acf(modelo3$residuals)
Ratifica-se o valor observado no teste de Durbin-Watson - vê-se
autocorrelações significantes
resid_1=rep(0,nEmp)
for(i in 2:nEmp)
resid_1[i]=modelo3$residuals[i-1]
modelo4= lm(tsEmp ~ t + t2 + factor(sazon) + resid_1)
summary(modelo4)
##
## Call:
## lm(formula = tsEmp ~ t + t2 + factor(sazon) + resid_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.52095 -0.13052 0.01346 0.14802 0.43917
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 90.4072543 0.1529018 591.276 < 2e-16 ***
## t -0.0471847 0.0070758 -6.668 1.98e-08 ***
## t2 0.0039730 0.0001044 38.042 < 2e-16 ***
## factor(sazon)2 0.1930471 0.1622457 1.190 0.240
## factor(sazon)3 0.1981482 0.1622749 1.221 0.228
## factor(sazon)4 0.7854122 0.1560452 5.033 6.63e-06 ***
## factor(sazon)5 2.0822974 0.1559728 13.350 < 2e-16 ***
## factor(sazon)6 3.5545698 0.1559378 22.795 < 2e-16 ***
## factor(sazon)7 3.8188962 0.1559400 24.490 < 2e-16 ***
## factor(sazon)8 3.4919433 0.1559799 22.387 < 2e-16 ***
## factor(sazon)9 2.9559550 0.1642073 18.001 < 2e-16 ***
## factor(sazon)10 1.9531826 0.1623179 12.033 2.22e-16 ***
## factor(sazon)11 1.5500678 0.1622725 9.552 7.34e-13 ***
## factor(sazon)12 0.1390069 0.1622454 0.857 0.396
## resid_1 1.0232201 0.0700189 14.613 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2565 on 50 degrees of freedom
## Multiple R-squared: 0.9974, Adjusted R-squared: 0.9967
## F-statistic: 1394 on 14 and 50 DF, p-value: < 2.2e-16
AIC(modelo4)
## [1] 22.53605
plot(modelo4)
shapiro.test(modelo4$residuals)
##
## Shapiro-Wilk normality test
##
## data: modelo4$residuals
## W = 0.98257, p-value = 0.4912
Ha evidencias para rejeitar H0 (H0: os residuos sao normalmente distribuidos)
durbinWatsonTest(modelo4)
## lag Autocorrelation D-W Statistic p-value
## 1 0.3417219 1.282271 0.002
## Alternative hypothesis: rho != 0
Não ha auto-correlacao de ordem 1
acf(modelo4$residuals)
Ratifica-se o valor observado no teste de Durbin-Watson - ve-se
autocorrelacoes significantes
lnEmp = log(tsEmp)
modelo5= lm(lnEmp ~ t + t2 + factor(sazon))
summary(modelo5)
##
## Call:
## lm(formula = lnEmp ~ t + t2 + factor(sazon))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.011754 -0.003275 -0.001043 0.003426 0.013831
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.503e+00 3.301e-03 1364.326 < 2e-16 ***
## t -2.003e-04 1.518e-04 -1.319 0.1929
## t2 3.601e-05 2.232e-06 16.138 < 2e-16 ***
## factor(sazon)2 2.116e-03 3.505e-03 0.604 0.5488
## factor(sazon)3 2.232e-03 3.505e-03 0.637 0.5271
## factor(sazon)4 8.636e-03 3.371e-03 2.562 0.0134 *
## factor(sazon)5 2.209e-02 3.369e-03 6.555 2.74e-08 ***
## factor(sazon)6 3.711e-02 3.368e-03 11.017 4.29e-15 ***
## factor(sazon)7 3.968e-02 3.368e-03 11.779 3.63e-16 ***
## factor(sazon)8 3.621e-02 3.369e-03 10.749 1.04e-14 ***
## factor(sazon)9 2.727e-02 3.508e-03 7.773 3.30e-10 ***
## factor(sazon)10 2.040e-02 3.506e-03 5.818 3.96e-07 ***
## factor(sazon)11 1.620e-02 3.505e-03 4.621 2.63e-05 ***
## factor(sazon)12 1.421e-03 3.505e-03 0.405 0.6869
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.005541 on 51 degrees of freedom
## Multiple R-squared: 0.9883, Adjusted R-squared: 0.9853
## F-statistic: 331.4 on 13 and 51 DF, p-value: < 2.2e-16
AIC(modelo5)
## [1] -476.7235
plot(modelo5)
shapiro.test(modelo5$residuals)
##
## Shapiro-Wilk normality test
##
## data: modelo5$residuals
## W = 0.98798, p-value = 0.7831
Ha evidencias para rejeitar H0 (H0: os residuos sao normalmente distribuidos)
durbinWatsonTest(modelo5)
## lag Autocorrelation D-W Statistic p-value
## 1 0.8105075 0.2567286 0
## Alternative hypothesis: rho != 0
Há evidencias para rejeitar-se H0 - H0: nao ha auto-correlacao de ordem 1
acf(modelo5$residuals)
Ratifica-se o valor observado no teste de Durbin-Watson - ve-se
autocorrelacoes significantes
resid_1=rep(0,nEmp)
for(i in 2:nEmp)
resid_1[i]=modelo5$residuals[i-1]
modelo6= lm(lnEmp ~ t + t2 + factor(sazon) + resid_1)
summary(modelo6)
##
## Call:
## lm(formula = lnEmp ~ t + t2 + factor(sazon) + resid_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.005222 -0.001443 0.000167 0.001492 0.004962
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.504e+00 1.537e-03 2930.595 < 2e-16 ***
## t -2.929e-04 7.098e-05 -4.127 0.000139 ***
## t2 3.774e-05 1.046e-06 36.071 < 2e-16 ***
## factor(sazon)2 2.089e-03 1.631e-03 1.281 0.206195
## factor(sazon)3 2.176e-03 1.631e-03 1.334 0.188363
## factor(sazon)4 8.466e-03 1.569e-03 5.397 1.86e-06 ***
## factor(sazon)5 2.190e-02 1.568e-03 13.967 < 2e-16 ***
## factor(sazon)6 3.690e-02 1.568e-03 23.540 < 2e-16 ***
## factor(sazon)7 3.945e-02 1.568e-03 25.164 < 2e-16 ***
## factor(sazon)8 3.596e-02 1.568e-03 22.933 < 2e-16 ***
## factor(sazon)9 3.002e-02 1.645e-03 18.252 < 2e-16 ***
## factor(sazon)10 2.046e-02 1.632e-03 12.536 < 2e-16 ***
## factor(sazon)11 1.624e-02 1.631e-03 9.956 1.87e-13 ***
## factor(sazon)12 1.444e-03 1.631e-03 0.885 0.380261
## resid_1 9.718e-01 7.136e-02 13.619 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002579 on 50 degrees of freedom
## Multiple R-squared: 0.9975, Adjusted R-squared: 0.9968
## F-statistic: 1434 on 14 and 50 DF, p-value: < 2.2e-16
AIC(modelo6)
## [1] -575.4465
plot(modelo6)
shapiro.test(modelo6$residuals)
##
## Shapiro-Wilk normality test
##
## data: modelo6$residuals
## W = 0.98736, p-value = 0.7498
Não ha evidencias para rejeitar H0 (H0: os residuos sao normalmente distribuidos)
durbinWatsonTest(modelo6)
## lag Autocorrelation D-W Statistic p-value
## 1 0.3358639 1.279857 0.004
## Alternative hypothesis: rho != 0
Há evidencias para rejeitar-se H0 - H0: nao ha auto-correlacao de ordem 1
acf(modelo6$residuals)
Ratifica-se o valor observado no teste de Durbin-Watson - ve-se
autocorrelacoes significantes
real=lnEmp[57:65]
n6=length(lnEmp)-9
serie=rep(0,n6)
for(i in 1:n6)
serie[i]=lnEmp[i]
serie=ts(serie,start=c(1999,04),frequency=12)
t=c(1:n6)
t2=t^2
sazon1=rep.int(4:12,1)
sazon2=rep.int(1:12,3)
sazon3=rep.int(1:11,1)
sazon6 = append(sazon1, sazon2)
sazon6 = append(sazon6, sazon3)
M5=lm(serie ~ t + t2+ factor(sazon6))
resid_m5=rep(0,n6)
for(i in 2:n6)
resid_m5[i]=M5$res[i-1]
M6= lm(serie ~ t + t2 + factor(sazon6) + resid_m5)
summary(M6)
##
## Call:
## lm(formula = serie ~ t + t2 + factor(sazon6) + resid_m5)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0048750 -0.0012220 -0.0001147 0.0011958 0.0052264
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.502e+00 1.629e-03 2763.130 < 2e-16 ***
## t 3.143e-05 8.442e-05 0.372 0.712
## t2 3.148e-05 1.437e-06 21.909 < 2e-16 ***
## factor(sazon6)2 2.185e-03 1.734e-03 1.260 0.215
## factor(sazon6)3 2.153e-03 1.734e-03 1.242 0.221
## factor(sazon6)4 8.771e-03 1.655e-03 5.300 4.25e-06 ***
## factor(sazon6)5 2.151e-02 1.653e-03 13.013 3.80e-16 ***
## factor(sazon6)6 3.607e-02 1.652e-03 21.829 < 2e-16 ***
## factor(sazon6)7 3.763e-02 1.652e-03 22.779 < 2e-16 ***
## factor(sazon6)8 3.310e-02 1.652e-03 20.038 < 2e-16 ***
## factor(sazon6)9 2.695e-02 1.653e-03 16.310 < 2e-16 ***
## factor(sazon6)10 2.013e-02 1.654e-03 12.173 3.38e-15 ***
## factor(sazon6)11 1.599e-02 1.656e-03 9.656 4.06e-12 ***
## factor(sazon6)12 1.968e-03 1.734e-03 1.135 0.263
## resid_m5 8.793e-01 7.456e-02 11.794 9.34e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002451 on 41 degrees of freedom
## Multiple R-squared: 0.9962, Adjusted R-squared: 0.9949
## F-statistic: 765.8 on 14 and 41 DF, p-value: < 2.2e-16
previsao = rep(1,9)
previsao[1] = M6$coef[1] + M6$coef[2]*(n6+1) + M6$coef[3]*(n6+1)^2 + M6$coef[14] + M6$coef[15]*M5$res[n6]
previsao[2] = M6$coef[1] + M6$coef[2]*(n6+2) + M6$coef[3]*(n6+2)^2 + M6$coef[15]*M5$res[n6]
for(i in 3:9)
previsao[i] = M6$coef[1] + M6$coef[2]*(n6+i) + M6$coef[3]*(n6+i)^2 + M6$coef[i+1] + (M6$coef[15])^i*M6$res[n6]
previsao
## [1] 4.608653 4.610338 4.615821 4.619523 4.629942 4.646553 4.665044 4.670609
## [9] 4.670153
exp(previsao)
## [1] 100.3489 100.5181 101.0708 101.4456 102.5082 104.2252 106.1703 106.7627
## [9] 106.7141
previsao=ts(previsao,start = c(2003,12), frequency = 12)
plot(exp(previsao))
real=ts(real,start = c(2003,12), frequency = 12)
plot(exp(real))
previsao
## Jan Feb Mar Apr May Jun Jul Aug
## 2003
## 2004 4.610338 4.615821 4.619523 4.629942 4.646553 4.665044 4.670609 4.670153
## Sep Oct Nov Dec
## 2003 4.608653
## 2004
real
## Jan Feb Mar Apr May Jun Jul Aug
## 2003
## 2004 4.609162 4.615121 4.620059 4.630838 4.652054 4.673763 4.685828 4.692265
## Sep Oct Nov Dec
## 2003 4.605170
## 2004
li=rep(1,9)
ls=rep(1,9)
previsaofinal=rep(1,9)
p=length(M6$coef)
PHI=1
for(i in 1:9)
{
li[i]=previsao[i] - qt(0.975,n6-p)*sd(M6$res)*sqrt(PHI)
ls[i]=previsao[i] + qt(0.975,n6-p)*sd(M6$res)*sqrt(PHI)
PHI=PHI+(M6$coef[15])^2*i
}
li
## [1] 4.604379 4.604646 4.608033 4.609372 4.617311 4.631382 4.647297 4.650265
## [9] 4.647200
for(i in 1:9)
{
previsaofinal[i] = exp(previsao[i])
li[i]=exp(li[i])
ls[i]=exp(ls[i])
}
length(previsaofinal)
## [1] 9
length(real)
## [1] 9
vetor_prev = cbind(exp(real),previsaofinal, li, ls)
vetor_prev
## exp(real) previsaofinal li ls
## Dec 2003 100.0 100.3489 99.92087 100.7788
## Jan 2004 100.4 100.5181 99.94755 101.0919
## Feb 2004 101.0 101.0708 100.28673 101.8610
## Mar 2004 101.5 101.4456 100.42111 102.4806
## Apr 2004 102.6 102.5082 101.22151 103.8112
## May 2004 104.8 104.2252 102.65580 105.8185
## Jun 2004 107.1 106.1703 104.30270 108.0713
## Jul 2004 108.4 106.7627 104.61275 108.9568
## Aug 2004 109.1 106.7141 104.29251 109.1918
for(i in 1:9)
SQP=(exp(real[i]) - previsaofinal[i])^2/9
SQP
## [1] 0.6325209
pp.test(lnEmp)
##
## Phillips-Perron Unit Root Test
##
## data: lnEmp
## Dickey-Fuller Z(alpha) = -9.6712, Truncation lag parameter = 3, p-value
## = 0.5413
## alternative hypothesis: stationary
Nao tem evidencias para rejeitar a hipotese nula H0 (H0: a série tem raiz unitária, logo é estacionária)
adf.test(lnEmp)
##
## Augmented Dickey-Fuller Test
##
## data: lnEmp
## Dickey-Fuller = -3.626, Lag order = 3, p-value = 0.03811
## alternative hypothesis: stationary
Rejeita-se a hipótese nula H0 (H0: a série não é estacionaria)
seguindo o teste Aumentado de Dickey-Fuller, A serie é estacionaria, entao necessita da componente Diferencial do ARIMA.
pacf(lnEmp)
observando o pacf, observa-se que ultrapassa-se a banda em ordem 1 - ou
seja, sera necessario um componente auto-regressivo 1
Box.test(lnEmp,lag=1,type="Ljung-Box")
##
## Box-Ljung test
##
## data: lnEmp
## X-squared = 56.314, df = 1, p-value = 6.173e-14
Box.test(lnEmp,lag=6,type="Ljung-Box")
##
## Box-Ljung test
##
## data: lnEmp
## X-squared = 205.6, df = 6, p-value < 2.2e-16
Box.test(lnEmp,lag=12,type="Ljung-Box")
##
## Box-Ljung test
##
## data: lnEmp
## X-squared = 317.18, df = 12, p-value < 2.2e-16
Box.test(lnEmp,lag=24,type="Ljung-Box")
##
## Box-Ljung test
##
## data: lnEmp
## X-squared = 354.4, df = 24, p-value < 2.2e-16
Box.test(lnEmp,lag=36,type="Ljung-Box")
##
## Box-Ljung test
##
## data: lnEmp
## X-squared = 423.84, df = 36, p-value < 2.2e-16
fit1 <- sarima(lnEmp,1,1,1,1,1,1,6)
## initial value -4.129883
## iter 2 value -5.283240
## iter 3 value -5.544635
## iter 4 value -5.563786
## iter 5 value -5.568725
## iter 6 value -5.578846
## iter 7 value -5.589521
## iter 8 value -5.590777
## iter 9 value -5.595660
## iter 10 value -5.596667
## iter 11 value -5.596768
## iter 12 value -5.596898
## iter 13 value -5.597348
## iter 14 value -5.597612
## iter 15 value -5.597809
## iter 16 value -5.598290
## iter 17 value -5.598463
## iter 18 value -5.599008
## iter 19 value -5.599301
## iter 20 value -5.599320
## iter 21 value -5.599350
## iter 22 value -5.599390
## iter 23 value -5.599402
## iter 24 value -5.599403
## iter 25 value -5.599403
## iter 26 value -5.599404
## iter 26 value -5.599404
## iter 26 value -5.599404
## final value -5.599404
## converged
## initial value -5.488519
## iter 2 value -5.490619
## iter 3 value -5.493312
## iter 4 value -5.496938
## iter 5 value -5.497854
## iter 6 value -5.498058
## iter 7 value -5.498195
## iter 8 value -5.498292
## iter 9 value -5.498507
## iter 10 value -5.498595
## iter 11 value -5.498616
## iter 12 value -5.498618
## iter 13 value -5.498618
## iter 14 value -5.498619
## iter 15 value -5.498619
## iter 16 value -5.498620
## iter 17 value -5.498620
## iter 18 value -5.498620
## iter 18 value -5.498620
## iter 18 value -5.498620
## final value -5.498620
## converged
## <><><><><><><><><><><><><><>
##
## Coefficients:
## Estimate SE t.value p.value
## ar1 0.7912 0.1807 4.3797 0.0001
## ma1 -0.2890 0.2567 -1.1258 0.2652
## sar1 -0.9040 0.0516 -17.5229 0.0000
## sma1 -0.5544 0.2138 -2.5935 0.0122
##
## sigma^2 estimated as 1.24632e-05 on 54 degrees of freedom
##
## AIC = -7.986949 AICc = -7.973936 BIC = -7.809324
##
fit1
## $fit
##
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S),
## include.mean = !no.constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ma1 sar1 sma1
## 0.7912 -0.2890 -0.9040 -0.5544
## s.e. 0.1807 0.2567 0.0516 0.2138
##
## sigma^2 estimated as 1.246e-05: log likelihood = 236.62, aic = -463.24
##
## $degrees_of_freedom
## [1] 54
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.7912 0.1807 4.3797 0.0001
## ma1 -0.2890 0.2567 -1.1258 0.2652
## sar1 -0.9040 0.0516 -17.5229 0.0000
## sma1 -0.5544 0.2138 -2.5935 0.0122
##
## $ICs
## AIC AICc BIC
## -7.986949 -7.973936 -7.809324
ma1 nao representativo, demais componentes respeitam as condicoes e sao representativos
Mantendo a mesma condicao de sazonalidade e removendo o componente media movel 1 (ma1) do modelo - indicado pelos resultados acima
fit2 <- sarima(lnEmp,1,1,0,1,1,1,6)
## initial value -4.129883
## iter 2 value -5.469381
## iter 3 value -5.566620
## iter 4 value -5.571658
## iter 5 value -5.576403
## iter 6 value -5.578365
## iter 7 value -5.578528
## iter 8 value -5.578827
## iter 9 value -5.578846
## iter 10 value -5.578851
## iter 11 value -5.578854
## iter 12 value -5.578859
## iter 12 value -5.578859
## iter 12 value -5.578859
## final value -5.578859
## converged
## initial value -5.475502
## iter 2 value -5.478572
## iter 3 value -5.482324
## iter 4 value -5.488341
## iter 5 value -5.489446
## iter 6 value -5.489628
## iter 7 value -5.489638
## iter 8 value -5.489640
## iter 8 value -5.489640
## final value -5.489640
## converged
## <><><><><><><><><><><><><><>
##
## Coefficients:
## Estimate SE t.value p.value
## ar1 0.5956 0.1249 4.7683 0.0000
## sar1 -0.8863 0.0539 -16.4451 0.0000
## sma1 -0.5848 0.2012 -2.9066 0.0053
##
## sigma^2 estimated as 1.27053e-05 on 55 degrees of freedom
##
## AIC = -8.003472 AICc = -7.995809 BIC = -7.861373
##
fit2
## $fit
##
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S),
## include.mean = !no.constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 sar1 sma1
## 0.5956 -0.8863 -0.5848
## s.e. 0.1249 0.0539 0.2012
##
## sigma^2 estimated as 1.271e-05: log likelihood = 236.1, aic = -464.2
##
## $degrees_of_freedom
## [1] 55
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.5956 0.1249 4.7683 0.0000
## sar1 -0.8863 0.0539 -16.4451 0.0000
## sma1 -0.5848 0.2012 -2.9066 0.0053
##
## $ICs
## AIC AICc BIC
## -8.003472 -7.995809 -7.861373
shapiro.test(fit2$fit$residuals)
##
## Shapiro-Wilk normality test
##
## data: fit2$fit$residuals
## W = 0.96492, p-value = 0.06214
Não ha evidencias para rejeitar a hipotese nula de que os residuos sao normalmente distribuidos
Observa-se melhora no AIC, aumento dos graus de liberdade e menor complexidade do modelo - este deve ser usado para simulacao
serie
## Jan Feb Mar Apr May Jun Jul Aug
## 1999 4.511958 4.526127 4.539030 4.541165 4.536891
## 2000 4.503137 4.506454 4.506454 4.514151 4.528289 4.548600 4.555980 4.554929
## 2001 4.518522 4.521789 4.522875 4.530447 4.541165 4.555980 4.558079 4.553877
## 2002 4.535820 4.541165 4.546481 4.557030 4.575741 4.591071 4.594109 4.593098
## 2003 4.574711 4.578826 4.579852 4.587006 4.601162 4.619073 4.621044 4.618086
## Sep Oct Nov Dec
## 1999 4.526127 4.519612 4.517431 4.503137
## 2000 4.549657 4.537961 4.534748 4.520701
## 2001 4.549657 4.548600 4.544358 4.533674
## 2002 4.591071 4.585987 4.584967 4.574711
## 2003 4.619073 4.619073 4.619073
serieSarima <- sarima.for(serie, 9, 1, 1, 0, 1, 1, 1, 6)
serieSarima$pred
## Jan Feb Mar Apr May Jun Jul Aug
## 2003
## 2004 4.609670 4.614122 4.613904 4.618848 4.631834 4.647129 4.649037 4.646800
## Sep Oct Nov Dec
## 2003 4.608695
## 2004
exp(serieSarima$pred)
## Jan Feb Mar Apr May Jun Jul Aug
## 2003
## 2004 100.4510 100.8992 100.8772 101.3772 102.7023 104.2851 104.4843 104.2508
## Sep Oct Nov Dec
## 2003 100.3531
## 2004
###Forecasting MSE
for(i in 1:9)
SQP=(exp(real[i]) - exp(serieSarima$pred[i]))^2/9
SQP
## [1] 2.61274
Dado tudo que foi visto ate o momento, pode-se afirmar que a serie contem uma componente de Nivel, de Tendencia e Sazonalidade, logo
fitHT <- HoltWinters(lnEmp)
plot(fitted(fitHT))
plot(fitHT)
fitHT2 <- HoltWinters(serie)
plot(fitted(fitHT2))
plot(fitHT2)
predHT2 = predict(fitHT2,9, prediction.interval = TRUE)
exp(predHT2)
## fit upr lwr
## Dec 2003 99.83901 100.6652 99.01956
## Jan 2004 99.77359 100.7571 98.79969
## Feb 2004 100.18276 101.3186 99.05965
## Mar 2004 100.38318 101.6638 99.11871
## Apr 2004 101.26748 102.6992 99.85568
## May 2004 102.81448 104.4077 101.24562
## Jun 2004 104.60144 106.3628 102.86927
## Jul 2004 105.04856 106.9577 103.17351
## Aug 2004 105.10032 107.1504 103.08943
exp(real)
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2003 100.0
## 2004 100.4 101.0 101.5 102.6 104.8 107.1 108.4 109.1
plot(exp(predHT2))
plot(exp(real))
for(i in 1:9)
SQP=(exp(real[i]) - exp(predHT2[i,1]))^2/9
SQP
## fit
## 1.77749
lnEmp
## Jan Feb Mar Apr May Jun Jul Aug
## 1999 4.511958 4.526127 4.539030 4.541165 4.536891
## 2000 4.503137 4.506454 4.506454 4.514151 4.528289 4.548600 4.555980 4.554929
## 2001 4.518522 4.521789 4.522875 4.530447 4.541165 4.555980 4.558079 4.553877
## 2002 4.535820 4.541165 4.546481 4.557030 4.575741 4.591071 4.594109 4.593098
## 2003 4.574711 4.578826 4.579852 4.587006 4.601162 4.619073 4.621044 4.618086
## 2004 4.609162 4.615121 4.620059 4.630838 4.652054 4.673763 4.685828 4.692265
## Sep Oct Nov Dec
## 1999 4.526127 4.519612 4.517431 4.503137
## 2000 4.549657 4.537961 4.534748 4.520701
## 2001 4.549657 4.548600 4.544358 4.533674
## 2002 4.591071 4.585987 4.584967 4.574711
## 2003 4.619073 4.619073 4.619073 4.605170
## 2004
model<-function(u){
mod<-dlmModSeas(frequency=12,dV=0,dW=c(exp(u[4]),rep(0,10)))+
dlmModPoly(2,dV=exp(u[1]),dW=(exp(u[2:3])))
}
outmle=dlmMLE(tsEmp,parm=rep(0,4),model)
exp(outmle$par)
## [1] 1.667719e-06 3.801566e-02 1.067618e-02 3.517537e-08
mod=model(outmle$par)
outmodFil=dlmFilter(tsEmp,mod)
outF<-dlmFilter(tsEmp,mod)
outS<-dlmSmooth(tsEmp,mod)
par(mfrow=c(4,1))
ts.plot((outF$m[10:nEmp,1]))
title("Sazonality: Filtering")
ts.plot((outF$m[10:nEmp,12]))
title("Slope: Filtering")
ts.plot((outF$m[10:nEmp,13]))
title("Level: Filtering")
ts.plot((outF$f[10:nEmp]),ylab="Yt^",xlab="t")
title("Preditos")
par(mfrow=c(1,1))
myt=matrix(NA,nEmp,2)
myt[,1]=tsEmp
myt[,2]=outF$f[1:nEmp]
ts.plot(myt,ylab="Yt^",xlab="t",col=c("black","red"))
Observa-se um descasamento nos preditos em momentos iniciais, depois o
modelo ajusta para as observações
ts.plot((outS$s[10:nEmp,1]))
title("Sazonality: Smoothing")
ts.plot((outS$s[10:nEmp,12]))
title("Slope: Smoothing")
ts.plot((outS$s[10:nEmp,13]))
title("Level: Smoothing")
# Análise dos resíduos
par(mfrow=c(1,1))
qqnorm(residuals(outF,sd=FALSE))
qqline(residuals(outF,sd=FALSE))
tsdiag(outF)
shapiro.test(residuals(outF,sd=FALSE))
##
## Shapiro-Wilk normality test
##
## data: residuals(outF, sd = FALSE)
## W = 0.94949, p-value = 0.009997
Ha evidencias para rejeitar H0 (H0: os residuos do modelo são normalmente distribuidos)
seriePura = exp(serie)
seriePura
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1999 91.1 92.4 93.6 93.8 93.4 92.4 91.8 91.6 90.3
## 2000 90.3 90.6 90.6 91.3 92.6 94.5 95.2 95.1 94.6 93.5 93.2 91.9
## 2001 91.7 92.0 92.1 92.8 93.8 95.2 95.4 95.0 94.6 94.5 94.1 93.1
## 2002 93.3 93.8 94.3 95.3 97.1 98.6 98.9 98.8 98.6 98.1 98.0 97.0
## 2003 97.0 97.4 97.5 98.2 99.6 101.4 101.6 101.3 101.4 101.4 101.4
outmle2=dlmMLE(seriePura,parm=rep(0,4),model)
exp(outmle2$par)
## [1] 1.948774e-06 4.435990e-02 3.268541e-03 1.822108e-07
mod2=model(outmle2$par)
outmodFil2=dlmFilter(seriePura,mod2)
outF2<-dlmFilter(seriePura,mod2)
outS2<-dlmSmooth(seriePura,mod2)
prev2=dlmForecast(outmodFil2,n=9)
prev2$f
## Jan Feb Mar Apr May Jun Jul Aug
## 2003
## 2004 100.7009 101.2891 101.6691 102.6409 104.2304 106.0133 106.5496 106.4995
## Sep Oct Nov Dec
## 2003 100.4796
## 2004
previsaoMEB=as.numeric(prev2$f)
lsMEB=previsaoMEB+1.96*sqrt(as.numeric(prev2$Q[1:9]))
liMEB=previsaoMEB-1.96*sqrt(as.numeric(prev2$Q[1:9]))
vetorPrevMEB = cbind(exp(real),previsaoMEB, liMEB, lsMEB)
for(i in 1:9)
SQP=(exp(real[i]) - previsaoMEB[i])^2/9
SQP
## [1] 0.7514231
A previsão obteve bons resultados com a Media do erro quadratico menor que 1